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Abstract 

The relative importance of Lagrangian and population dynamics on spatial pattern 
formation in the distribution of plankton near the ocean's surface is investigated. 
Phytoplankton and zooplankton are treated as biologically interacting tracers that 
are passively advected by a divergence-free two-dimensional nonsteady velocity field. 
The time-dependence of this field is assumed to be quasiperiodic, which extends the 
analysis presented beyond standard chaotic advection. Various forms of predator- 
prey interactions, including complex interactions, are considered. Numerical sim- 
ulations illustrate how Lagrangian dynamics, which is characterized by a mixed 
phase space with coexisting regular and chaotic motion regions, influences plankton 
patchiness generation depending on the nature of the biological interactions. 
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1 Introduction 



The central role played by the peculiar topology of the phase space associ- 
ated with the Lagrangian dynamics of a deterministic divergence-free two- 
dimensional oscillatory velocity field in controlling dispersion of passive trac- 
ers on the surface of the ocean was first emphasized by Pasmanter In that 
work a kinematic model for advection by a tidal current is used to provide a 
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possible explanation for the observed [2j anomalous (i.e., non-Brownian) dif- 
fusion in diagrams of the average square size of pollutant clouds as a function 
of time in shallow seas. Because oscillations of the velocity field in 0] are as- 
sumed to be periodic, the Lagrangian motion is governed by a one-and-a-half- 
degree-of-freedom Hamiltonian system. The phase space of this motion, whose 
coordinates are those of the physical space for advection, is characterized by 
the coexistence of regular and chaotic fluid particle trajectories. More specifi- 
cally, the phase space exhibits a mixed structure composed essentially of three 
types of topological objects: chaotic motion regions; invariant tori — so-called 
KAM tori after Kolmogorov, Arnold and Moser — which carry quasiperiodic 
motion and separate the chaotic regions; and other invariant tori — so-called 
resonance islands — which are built up also from quasiperiodic trajectories and 
populate the chaotic regions [cf., e.g.,0]. The presence of the invariant sets of 
regular motion fundamentally influences the dispersion properties of clouds of 
points (fluid particles) in phase space [cf., e.g.,|4j by two mechanisms. First, 
both KAM tori and resonance islands serve as impenetrable barriers for trans- 
port. Second, phase-space points can experience stickiness on the boundary 
of the resonance islands, and also undergo long excursions (so-called Levy 
flights) within the chaotic regions. As a consequence, the dispersion of clouds 
of phase-space points deviates from a normal (i.e., Brownian) diffusive pro- 
cess, in accordance with the observed anomalous dispersive characteristics of 
passive tracers. 

An important observation is that the above type of Lagrangian dynamics, 
which we refer to as Lagrangian mixed phase-space dynamics, is not exclusively 
restricted to advection by a time-periodic velocity field, which corresponds to 
the standard chaotic advection paradigm . Rather, velocity fields with com- 
plicated time (and also spatial) dependence have been shown 0, 0, to pro- 
duce Lagrangian mixed phase-space dynamics as well. By complicated time 
dependence we specifically mean quasiperiodic time dependence where the 
number of basic frequencies is very large. We remark, however, that the above 
observation should not come as a surprise inasmuch as KAM theory [cf., e.g., 
[if has been extended 0, 0] to include near-integrable nonautonomous Hamil- 
tonians with quasiperiodic time dependence. Such time dependence arises nat- 
urally when a highly structured velocity field is represented in Fourier space 
as the superposition of many harmonics, which may include randomly cho- 
sen phases, and whose amplitudes are taken from an appropriate wavenumber 
spectrum with some appropriately chosen low- and high-wavenumber cutoffs. 
This indicates, for instance, that small-scale environmental perturbations to a 
large-scale velocity field can be treated deterministically without fundamen- 
tally changing the Lagrangian mixed phase-space dynamics resulting from 
consideration of the large-scale advection component in isolation. Further- 
more, this suggests a quasi-deterministic representation of mesoscale turbu- 
lence, thereby leading to the expectation that Lagrangian mixed phase-space 
dynamics may also play a role in open-ocean environments. The latter expec- 
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tation finds some support in the analysis of surface drifter and submerged float 
trajectories [e.g., ll[ H, 3, li> HI; which have revealed fractal properties. 

From the perspective of the distribution of a passive tracer (i.e., any scalar that 
does not alter the advection field) the chaotic motion regions of a mixed phase 
space correspond to regions of vigorously stirred tracer. The regular motion 
sets, particularly the resonance islands, correspond to patches of poorly stirred 
tracer. For biologically interacting passive tracers, such as phytoplankton and 
zooplankton, the relationship between the spatial patterns observed in their 
distributions on the ocean's surface [e.g., [l(J and the features that comprise 
a mixed phase space is not obvious. For instance, the effect of fluid particle 
clustering in a resonance island may be compensated by the occurrence of 
nonequilibrated population dynamics leading to a region of heterogeneously 
distributed plankton rather than a patch of homogeneous plankton concentra- 
tion. Numerical evidence [e.g.,[l7, 3, 1^, 2(1 U, 22] nonetheless suggests that 



advection by ocean currents can indeed play an important role in controlling 
the generation of spatial patterns in plankton distributions. However, the nu- 
merical evidence presented have typically relied on either the assumption of 
equilibrated population dynamics or the consideration of a specific biological 
interaction scenario. So far as we are aware, no work has been reported on 
plankton patchiness generation that focuses on the competing roles played 
by Lagrangian mixed phase-space dynamics and the predator-prey dynamics 
controlling the plankton density intrinsic evolution. 

The main goal of the present contribution is to carry out such a study. To keep 
our analysis simple, we will consider a kinematic model of advection bv a tidal 
current that is only somewhat more complicated than that used in More 
precisely, we will take into account the spring-to-neap cycle associated with 
the superposition of the semilunar (semidiurnal lunar) and semisolar (semidi- 
urnal solar) tides, which leads to a streamfunction (Hamiltonian) that depends 
quasiperiodically on time. Consideration of this time dependence will extend 
the scope of our work beyond the realm of standard chaotic advection, which 
differs in some respects from advection by a time-quasiperiodic velocity field. 
A secondary goal of the present work is discuss these differences, which do not 
depend on the number n of basic frequencies considered, thereby allowing us 
to choose n = 2 with no loss of generality. Unlike prior studies of plankton 
patchiness generation, we will consider a variety of biological interaction types, 
including complex interactions, which will be assumed to be governed by a 
simple set of predator-prey equations. Specifically, in addition to autonomous 
predator-prey interactions we will also consider nonautonomous predator- 
prey dynamics by allowing for a diurnal cycle in the food availability to prey, 
which will lead to a very rich set of interaction scenarios, including chaotic 
predator-prey interactions. Remarkably, even though the biological relevance 
of the latter interactions is supported both theoretically and experimentally 
[jjj and references therein], they have been excluded in prior studies of the 
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production of plankton patchiness. The rationale for their exclusion may be 
that their presence does not lead to the development of spatial patterns in the 
plankton distributions. However, anticipating some of our results, in this paper 
we will show that Lagrangian mixed phase-space dynamics can lead to plank- 
ton patchiness development even under chaotic predator-prey interactions. It 
will be shown that in general the influence of Lagrangian mixed phase-space 
dynamics plays on plankton patchiness generation is less constrained by the 
type of biological interactions during early stages of the evolution than during 
later stages. 

The remainder of the paper has the following organization. In Section 2 we 
introduce the time-quasiperiodic kinematic velocity field model that is used 
throughout the paper. In the same section we discuss the effects of the time- 
quasiperiodicity of the velocity field on Lagrangian dynamics. This is done 
by suitably constructing a Poincare section. Section 3 discusses the dynamics 
of tracers passively advected by the velocity field. Section 4 is devoted to 
study the dynamics of phytoplankton and zooplankton by treating them as 
biologically interacting tracers that are passively advected by the velocity field 
defined earlier. The different types of biological interactions considered, both 
autonomous and nonautonomous, are discussed prior to discussing the effects 
of advection. Section 5 offers a discussion on several topics including the effects 
of random fluctuations of the velocity field and the relevance of the present 
work to the study of plankton advection by more complicated ocean currents. 
A summary and the conclusions of the paper are finally presented in Section 
6. 



2 Fluid Particle Dynamics 



We consider a divergence-free velocity field u = (d y tp, —d x ip) defined at each 
time t on the Euclidean plane with Cartesian coordinates (x, y) by a stream- 
function ip of the form |l| 

ij;(x, y, t) := u(t)y - v(t)x + U(t)l~ x cos ly - V{t)k~ l cos kx, (2.1) 

where 

n 

u(t) = uo + Ui cos o"jt (2.2) 
i 

and similarly for v(t), U(t), and V(t). The velocity field u includes the essential 
features of a tidal (i.e., oscillatory) velocity field in a shallow sea wherein ebb 
and flood surpluses organize in cell-like circulation structures. A dynamical 
explanation of these "tidal residual eddies" can be found in 24 , |25| . In £j| 



the kinematic model (2.1) is written for arbitrary n with the frequencies {a{\ 
chosen to be commensurate, which leads to periodic time dependence. Analysis 
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in |l| is restricted to the n = 1 case, which corresponds to the standard 
chaotic advection paradigm. A model similar to (2.1) with n = 1 is also 
used in |26|. In this paper we consider a time-quasiperiodic dependence with 
n = 2 incommensurate frequencies. This choice presents all the mathematical 
difficulties associated with a large number of incommensurate frequencies, and 
furthermore extends our analysis beyond standard chaotic advection. 

More precisely, we choose o\ = 27r/12.4206 rad h _1 and 02 = rad 
h _1 , which presumes a primarily semidiurnal tidal regime that includes the 
spring-to-neap tidal cycle associated with the superposition of the semilunar 
tidal component and a smaller semisolar tidal component. For this frequency 
choice u is actually time-periodic with a period of exactly 1495224 h (rs 170 
y); however, this period is too long and u is more appropriately interpreted as 
being time-quasiperiodic. In |26| it is speculated about the possible similarities 
of the effects of this time-quasiperiodicity on fluid particle motion compared 
to the effects of time-periodicity. In this paper we show that there indeed are 
qualitative similarities but also quantitative differences. 

The trajectories (x(t),y(t)) of fluid particles obey 

x = d y if), y = -d x ip. (2.3) 

Upon identifying x—y with the canonically conjugate coordinate-momentum 
pair, system (2.3) is seen to constitute a one-degree-of-freedom Hamiltonian 
system that is nonautonomous because the Hamiltonian, ip, depends explicitly 
on t. More specifically, such dependence is assumed of the form 



if)(x,y,t) = *(x,y,a 1 t,a 2 t), (2.4) 

where * : E 2 x T 2 -> E 2 (T = E/2vrZ is the standard one-torus). The 
phase space of (2.3) is constituted by the real two-space E 2 endowed with the 
standard symplectic form dx A dy, which is preserved by the time-dependent 
flow of u. Namely, the family of (t, s)-advance mappings of E 2 to itself 



<f> ttS :{x{t),y(t))^{x{s),y(s)) (2.5) 

satisfying 4> ts o(fi s s+r = (p t s+r and <f> t t = id. Invariance of dx A dy under 
the flow of u is equivalent to preservation of area in phase space (Liouville's 
theorem) . 



The time-independent part of if) is integrable; typical fluid particle trajectories 
are depicted in the left panel of Fig. 1 on E/27r/c _1 Z x E/27r/ _1 Z. In general, 
due to explicit time dependence, if) is not integrable and trajectories may be 
chaotic, i.e., exhibit sensitivity to initial conditions. A two-dimensional por- 
trait of the dynamics can be sketched by constructing a second-order Poincare 
section [cf., e.g.,|2j], Chapter 2.3]. This is accomplished by taking the first re- 
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Fig. 1. Phase portraits of (2.3) with ip defined by (2.1) excluding time-dependent 
components (left) and including two time-dependent components (right) with fre- 
quencies, which, for any practical purpose, can be interpreted as being incommen- 
surate. Parameter choices are: a\ = 27r/12.4206 rad h _1 , 02 = 27r/12 rad h _1 , 
k = 1 = 2vr/10 km" 1 , u = -0.5 km d -1 , v = km d -1 , U = V = 1 km d -1 , 
u\ = v 1 = 112 = V2 = —0.1 km d _1 , and U\ = V\ = U2 = V2 = 0.1 km d _1 . The plot 
in the right panel corresponds to a second-order Poincare section of (2.3) defined 
by (2.6) with 5 = 10" 4 . 

turn map under the flow of (2.3) to the set 



where (f { := a it mod27r and 5 > is small, and then plotting the resulting 
points using their E 2 coordinates. 1 

The right panel of Fig. 1 shows a second-order Poincare section of (2.3) as- 
suming 5 = 10 -3 . To ensure strict Hamiltonian behavior, the flow of (2.3) 
was numerically computed using a fully-implicit sixth-order Gauss-Legendre 
method [cf., e.g.,|28|, which is symplectic (i.e., area preserving). Despite the 
blurring associated with finiteness of S, the right panel of Fig. 1 reveals quite 
clearly a mixed phase space with coexisting regular and chaotic motion re- 
gions. Such a mixed phase space is most commonly associated with Poincare 
sections of (2.3) with time-periodic ip constructed by strobing trajectories at 
the frequency of ip. Legitimacy of the mixed phase-space topology identified 
here with the aid of the second-order Poincare section technique is insured by 
the rigorous KAM-theory- related results of Jorba and Simo |9]. 

To grasp the meaning of the mixed phase-space dynamics associated with (2.3) 
for time-quasiper iodic ?/>, it is most convenient to consider ip split into back- 



Because 5 can never be chosen to be equal to zero, the set S actually lies on 
a three-dimensional space. Consequently, a second-order Poincare section is not a 
plane in a strict sense; rather, it is a slice of nonzero, albeit thin, thickness. 



E := R 2 x {ip x = 0, y? 2 e [0,5]} 



(2.6) 
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ground, if) (I), and perturbation, ■0 1 (J,i?,i) + ip 2 (I ,-d,t) where i/j^I, t) = 
i?, (Jit). Here I-d are some locally defined action-symplectic-angie vari- 
ables. When ip 1 + ip 2 — 0) motion lies on one-tori {/ = const} with frequency 
cu(I) := 4>' (I) (lighter curves in the left panel of Fig. 1). When + ijj 2 \ ^ 
is sufficiently small, the KAM theorem proved in |j| guarantees that the back- 
ground tori for which the frequencies {oo{I), o"i, cr 2 } are incommensurate in a 
Diophantine sense survive the perturbation, albeit slightly distorted and with 
the frequencies <7i and <r 2 added to the unperturbed frequency The latter 
result, which applies to almost all background tori, holds under the assump- 
tions that ip is sufficiently smooth, o\ and <t 2 satisfy a Diophantine inequality, 
and ijj is nondegenerate in the sense to' (I) ^ 0. 2 The main effect of the per- 
turbation is to impart on a background torus a quasiperiodic vibration with 
frequencies <Ti and <r 2 , which has been referred in jjjj to as a "quasiperiodic 
dancing." This contrasts with the widely studied time-periodic perturbation 
case for which the vibration is periodic. 

The open and concentric curves in the right panel of Fig. 1 are examples of 
KAM tori 3 of the above quasiperiodically vibrating type. The closed curves 
flanking the KAM tori correspond to tori into which background tori have bro- 
ken up due to resonance, e.g., u(Io)/ai =p/g e Q. These tori (so-called reso- 
nance islands) enclose elliptic points separated by an equal number of saddle 
points, which, unlike the well-known time-periodic perturbation case, corre- 
spond to quasiperiodic trajectories, e.g., with frequencies {qu(I ) = pai, cr 2 }. 
Associated with the saddle points are stable and unstable manifolds that trans- 
versely intersect infinitely many times producing a narrow band of chaotic 
motion (so-called chaotic layer) as do perturbed heteroclinic trajectories [e.g. 
33|. The fuzzy two-dimensional regions filled up with dots in the right panel of 



Fig. 1 correspond to chaotic motion regions. These regions originate from the 
merging of the chaotic layers associated with the breakup of two homoclinic 
background trajectories (darker curves in the left panel of Fig. 1) and neigh- 
boring chaotic layers associated with the breakup of resonant background tori 
that are not separated by KAM tori. 



The above results straightforwardly extend [9| to the case of an arbitrary 
number n of perturbation frequencies. When n is large, however, other phase- 
space visualization techniques must be considered. Examples are the estima- 
tion of Lyapunov exponents and a variant of the Poincare section technique 
designed in [34] to detect domains of finite-time stability in the phase space 



However, there are conditions [e.g.,0 on higher-order derivatives of io(I) that 
probably 30] allow to make density estimates of tori even in the neighborhood of a 
point where to' {I) = 0. 

3 More precisely, cantori, due to insufficient incommensurability of a\ and o"2- Can- 
tori are gappy tori (more rigorously, Cantor sets) on which the dynamics is similar 
to that on KAM t ori, except that they can be occasionally traversed by trajectories 
through the gaps 0, . 
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of randomly perturbed systems. In [6( two-dimensional charts of finite-time 
estimates of Lyapunov exponents were constructed for the Lagrangian motion 
produced by a gyre-like ocean circulation composed of a steady component 
and a unsteady perturbation with a broad band of frequencies. These charts 
revealed regions of small absolute numerical values, suggesting KAM tori and 
resonance islands, and regions of larger absolute numerical values, identifi- 
able with chaotic motion regions. Similar results were found 0, HI using the 
Poincare section technique of applied to the Lagrangian motion produced 
by an oceanic jet interacting with an eddy in a noisy environment. This numer- 
ical evidence suggests that the results presented below also apply to the large 
n case, which allows to represent more complicated (e.g., turbulent) advection 
fields (cf. Section 4). 



3 Passive Tracer Dynamics 

The concentration of a passive tracer 9(x,y,t) evolves according to the Liou- 
ville (advection) equation 

dt9 + [e,ip] = 0, (3.1) 

where [f,g] '■— d x fd y g — d y fd x g denotes the symplectic bracket with respect 
to da; A dy. Equation (3.1) is equivalent to 6 = along (x(t),y(t)) satisfying 
(2.3), i.e., x = [x,ip] and y = [y,if)]. Consequently, the solution of (3.1) for 
certain initial distribution 9(x,y, 0) may be expressed as 

e(x,y,t) = 9((f> Q) _ t (x,y),0). (3.2) 

In other words, the distribution of 9 at any instant is an area-preserving rear- 
rangement of its initial distribution. 

Existence and uniqueness of solutions of (2.3) guarantee that trajectories with 
different initial conditions can never intersect. KAM tori and resonance islands 
are invariant objects in the sense that they are built up from trajectories; 
hence, they can never be traversed by trajectories with different initial condi- 
tions. In fluid mechanics jargon, these objects are said to be material meaning 
that they are composed by the same fluid particles. The chaotic stirring of a 
passive tracer is thus constrained by KAM tori, which are barriers for tracer 
transport, and resonance islands, which constitute actual tracer patches. Both 
types of mixed space-space structures lead to spatial pattern formation in the 
continuous-time distribution of a passive tracer. 

Fig. 2 illustrates the above for an arbitrary tracer that is passively advected 
by the quasiperiodic tidal current defined by streamfunction (2.1) with pa- 
rameters as in Fig. 1. To construct Fig. 2 we advected 250000 fluid particles 
distributed initially over a regular grid defined on a doubly-periodic 10 km x 
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Fig. 2. Snapshots of the evolution of a tracer that is passively advected by the 
time-quasiperiodic tidal velocity field defined by streamfunction (2.1) with the same 
parameters used to construct the Poincare section in the right panel of Fig. 1. The 
color scale indicates tracer concentration or simply a fluid particle labeling. 

10 km domain of R 2 . The spatial distribution of the tracer at any stage of the 
evolution is obtained from Delaunay triangulation of fluid particle positions 
by triangle-based cubic interpolation onto the regular grid. The distribution 
of the tracer at t = is assumed to be a Gaussian. At t = 15 d the tracer 
distribution reveals vigorously stirred regions along with poorly stirred regions 
in the form of spiral-like vortices and filaments. These features resemble quite 
closely those embodying the mixed phase space shown in the right panel of 
Fig. 1. The mixed phase-space features are much more clearly visible in the 
tracer distribution at the later instant t = 60 d, which can be considered long 
enough for (constrained) chaotic stirring to develop substantially. Note, how- 
ever, that 60 d is much shorter than the integration time required to construct 
the second-order Poincare section in the right panel of Fig. 1. 

An important observation is that the tracer distribution features can never 
match those of the mixed phase space because, as explained above, the mixed 
phase space features "dance" with the "rhythm" of the quasiperiodic tide. 
This contrasts with the standard chaotic advection paradigm (time-periodic 
if}) for which the motion of the mixed phase space structures is periodic. In 
that case matching can be achieved by displaying the tracer distribution at 
some (sufficiently large) multiple of the time-period of if}. 



4 Biologically Interacting Passive Tracer Dynamics 

This section is devoted to show that Lagrangian mixed phase-space dynamics 
can play a distinguished role in setting the distribution of biologically inter- 
acting tracers depending on the nature of the interaction. We will assume that 
the biologically interacting tracers obey a coupled set of Liouville equations 
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of the form 

a t e + [o,ii>] = F {e,t), (4.i) 

where 9 is understood as a vector of biologically interacting tracers and Fg(6, t) 
represents the interaction. The solution of (4.1) may be expressed as 

9(x, y, t) = 6 ieat ((j) 0i _ t (x, y),t), (4.2) 

where 9 rest (x, y, t) is the solution of (4.1) with initial distribution 8 rest (x, y, 0) = 
9(x, y, 0) assuming that the fluid is at rest. In other words, the distribution of 
9 at time t is given by an area-preserving rearrangement of the distribution it 
would have attained at that same instant if the fluid were at rest. 

Our particular focus will be on the plankton distribution problem, which we 
will approach by considering the interaction between phytoplankton, P, and 
zooplankton, Z, in a two- level food chain. We will consider nonautonomous 
Ro^nzweig-MacArthur predator -prey interactions governed by 

F P (P, Z, t) = r(t) (l - P - (4.3a) 



K(t)J P s + P' 

P + P _ 



ecPZ 

F z (P,Z,t) = ^—--dZ, (4.3b) 



where 



r(t) = r (1 - ecosfit), K(t) = K (l - ecos fit). (4.4) 

In this model the prey density, P, grows logistically with growth rate r(t) > 
and carrying capacity K(t) > in the absence of predation, Z = 0. The preda- 
tor consumes the prey according to a Holling's type-II functional response, 
which presumes a maximal consumption rate c > and half-saturation den- 
sity P s > 0. The model further assumes that the predator's per capita death 
rate, d > 0, is independent of density, and that its numerical response is pro- 
portional to its functional response with efficiency e > 0. We will set 2n/Q = 1 
d so that (4.4) crudely accounts for the diurnal cycle of light availability. 

Other, more sophisticated and perhaps more realistic, two- or higher-level 
food chain models and multispecies systems have been proposed, which could 
of course be considered. However, model (4.3), while very simple, describes 
a very rich variety of biological interaction types, which is enough for the 
purposes of the present work. 



4-1 Biological Interactions 



The predator-prey interactions governed by (4.3), i.e., P = Fp(P,Z,t) and 
Z = Fz(P, Z,t) along (x(t),y(t)), constitute a dissipative dynamical systems 
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provided P e /K ^ and P e /P s 7^ 0, where 



Pe := (4.5) 
ec — a 

which is biologically most sensible. Consequently, the topology of phase space 
for predator-prey interactions will be fundamentally different than the topol- 
ogy of phase space for fluid particle motion governed by (2.3) and discussed 
in Section 2. Predator-prey interactions are described in the next paragraph 
for a more detailed description the reader is referred to Refs. 



text par 
0, Kl- 



in the autonomous case (e = 0), predator-prey interactions obeying (4.3) 
possess a unique positive equilibrium (P e , Z e ), where 

Z e :=^(l-^W + P e ), (4.6) 



c V K 

provided < P c /K < 1 and P e /P s > 0. This equilibrium is subject to 
a change of stability type through a Hopf bifurcation. For P e /K > 1/(2 + 
P s P e ) predator and prey densities undergo damped oscillations toward (P e , Z e ), 
while for P e / K < 1/(2 + P s / P e ) predator and prey densities experience time- 
asymptotic, self- sustained oscillations about (P e ,Z e ). 

In the nonautonomous case (e > 0), the equilibrium (P e , Z e ) is replaced by 
a small-amplitude periodic orbit that loses stability through discrete Hopf or 
Neimark-S acker bifurcations. Chaotic predator-prey interactions, which are 
characterized by the presence of a strange attractor, eventually give rise via 
either torus destruction (quasiperiodic route to chaos) or a cascade of period 
doublings (subharmonic route to chaos). The first kind of strange attractors 
results for small e when the autonomous predator-prey interactions lie on a 
limit cycle. Strange attractors of the second class develop for larger e, even 
when the autonomous predator-prey system does not cycle. 

Fig. 3 depicts phase portraits for autonomous interactions that decay toward 
a stable spiral equilibrium (left panel) and asymptotically cycle (right panel). 
Fig. 4 shows phase portraits for nonautonomous interactions taking place on 
a periodic 1 : 18 subharmonic attractor (left panel) and lying on a strange 
attractor attained via quasiperiodic route (right panel). The latter are shown 
on Poincare sections of the predator-prey system, which are computed by 
strobing orbits at the environmental variability frequency Q/27T. To construct 
Figs. 3-4 the predator-prey equations were numerically integrated using the 
Dormand-Prince 5(4) algorithm. The period of the limit cycle oscillation in the 
right panel of Fig. 3 is of about 9 d. Estimated Lyapunov exponents scaled by Q 
for the strange attractor in the right panel of Fig. 4 are {+0.0239, 0, —0.0539}. 
This gives a predictability horizon (inverse of the positive Lyapunov exponent) 
not exceeding 10 d for the chaotic interactions considered here. 
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Fig. 3. Phase portraits for phytoplankton, P, and zooplankton, Z, interacting ac- 
cording to autonomous predator-prey dynamics governed by (4.3), i.e., with e = 
in (4.4). The thin-solid curve in the right panel depicts a damped oscillation to- 
ward (P e ,Z e ), defined in Eqs. (4.5)-(4.6), for lnP e /K = -3 and lnP e /P s = -4. 
The thick-solid closed curve in the left panel illustrates a self-sustained oscillation 
about the equilibrium (P e ,Z e ) for In P c /Kq = —2.8 and In P e /P s = —1.6; orbits for 
two different initial conditions are shown in this panel as thin-solid curves, which 
collapse into the limit cycle represented by the thick-solid curve. In both panels 
ro/Q = 0.4475 and d/Q, = 0.2011. The blue dots in the left panel correspond to a 
discrete set of P—Z values in the range spanned by the distributions shown in Figs. 
5-6 at t = 0. The green (resp., red) dots in the left and right panels correspond 
to a discrete set of P—Z values in the range spanned by the distributions shown, 
respectively, in the top and mid-top row panels of Figs. 5-6 at t = 15 d (resp., 
t = 60 d). 

Other, intermediate, forms of interaction (not shown but discussed below) 
are nonautonomous interactions lying on a torus or quasiperiodic attractor, 
and interactions possessing multiple periodic or quasiperiodic attractors. The 
former are excited for sufficiently small e when the autonomous interactions 
take place on a limit cycle; the latter precede chaotic interactions. 

4-2 Biological Interactions and Passive Advection 

Figs. 5 and 6 depict, respectively, possible evolution scenarios for phytoplank- 
ton and zooplankton, which, while interacting through predator-prey dynam- 
ics obeying (4.3), are passively advected by the same quasiperiodic tidal cur- 
rent involved in the construction of Figs. 1-2. Predator-prey interaction pa- 
rameters in the top and mid-top row panels of Figs. 5 and 6 are, respectively, as 
in the left and right panels of Fig. 3. Interaction parameters in the mid-bottom 
and bottom row panels of Figs. 5 and 6 are, respectively, as in the left and right 
panels of Fig. 4. Figs. 5-6 were constructed by advecting 250000 fluid parti- 
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Fig. 4. As in Fig. 3, but for nonautonomous predator-prey dynamics governed by 
(4.3), i.e., with e ^ in (4.3b). The left and right panels correspond to Poincare 
sections showing, respectively, a periodic attractor (1 : 18 subharmonic) attained 
with e = 0.2 and a strange attractor attained with e = 0.6; in both panels 
In P e /K = -2.8, lnP /P s = -1.6, r /n = 0.4475, and d/fl = 0.2011. The green 
(resp., red) dots in the left and right panels correspond to a discrete set of P-Z val- 
ues in the range spanned by the distributions shown, respectively, in the mid-bottom 
and bottom row panels of Figs. 5-6 at t = 15 d (resp., t = 60 d). 

cles, each one carrying phytoplankton and zooplankton densities, distributed 
on regular grid denned on a doubly-periodic 10 km x 10 km domain of R 2 
as in Fig. 2. The spatial distributions of phytoplankton and zooplankton at 
any stage of the evolution are obtained by triangle-based cubic interpolation 
onto the regular grid based on Delaunay triangulation of the fluid particle 
positions. 

At t = the phytoplankton distribution (Fig. 5, left column panels) is assumed 
to be a Gaussian, while the zooplankton distribution (Fig. 6, left column 
panels) is assumed to be uniform. This corresponds to a sudden (nonuniform) 
availability of prey resources (phytoplankton blooming) over the entire domain 
occupied by the predator. The locus in phase space of the initial plankton 
distributions is indicated by blue dots in Fig. 3-4. 

Distributions of phytoplankton and zooplankton after 15 d (resp., 60 d) are 
depicted in the middle (resp., right) column panels of Figs. 5 and 6, respec- 
tively. The locus in phase space of these distributions at t — 15 d (resp., t = 60 
d) is indicated in Figs. 3-4 with green (resp., red) dots. Note that t — 15 d 
is still too early for the regular autonomous predator-prey interactions in the 
left panel of Fig. 3 to reach an asymptotically stable equilibrium, or for the 
regular nonautonomous predator-prey interactions in the left panel of Fig. 4 to 
settle on a periodic attracting set. Likewise, at t = 15 d chaotic predator-prey 
interactions in the right panel of Fig. 4 produce phytoplankton and zooplank- 
ton density values that cover a smaller portion of the corresponding strange 
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attractor than those attained at t = 60 d. On the other hand, approximately 
15 d must elapse for the autonomous predator-prey interactions in the left 
panel of Fig. 3 to reach a limit-cycle oscillation. 

4-2.1 The Early Evolution Stage 

At t — 15 d phytoplankton (Fig. 5, middle column panels) and zooplankton 
(Fig. 6, middle column panels) show structures in their distributions that are 
very similar to those in the distribution of a passive tracer (Fig. 2, middle 
panel), which are an early manifestation of Lagrangian mixed phase-space dy- 
namics (Fig. 1, right panel). Remarkably, although the estimated predictabil- 
ity horizon for the chaotic predator-prey interactions is shorter than 15 d, 
Lagrangian mixed phase-space dynamics shows up very clearly in the plank- 
ton distributions. The reason for this is that sensitivity to initial conditions 
do not manifest at the same time for all orbits emanating from the latter 
in the phase space for predator-prey interactions. A careful examination of 
the mid panel in the bottom row of Figs. 5-6, however, reveals small regions 
of very irregularly distributed plankton interspersed within more uniformly 
distributed plankton, which are not a result of chaotic stirring but rather 
of chaotic predator-prey interactions. Other regions of very irregularly dis- 
tributed plankton resulting from chaotic interactions may lie hidden within 
the chaotically stirred plankton regions. Clearly, for these regions to get ex- 
posed the timescale for substantial development of chaotic stirring must be 
much larger than the predictability horizon of chaotic predator-prey interac- 
tions. 



4-2.2 The Late Evolution Stage 

At t = 60 d the regular predator-prey interactions in the left panel of Fig. 
3 have almost reached an asymptotically stable equilibrium. Accordingly, the 
right panel in the top row of Figs. 5 and 6 show very uniformly distributed 
phytoplankton and zooplankton concentrations, respectively. 

At the same instant t = 60 d, regions of vigorously stirred phytoplankton 
and zooplankton concentrations along with less stirred regions in the form 
of filaments and spiral-like vortices develop when regular predator-prey in- 
teractions autonomously cycle (Figs. 5-6, right panel in the mid-top row) or 
nonautonomously behave on a 1 : 18 subharmonic attractor (Figs. 5-6, right 
panel in the mid-bottom row). These features resemble very closely those of 
the passive tracer distribution at t = 60 d in the right panel of Fig. 2, which 
are a manifestation of the mixed phase-space topology that characterizes La- 
grangian dynamics. However, note that for predator-prey interactions lying 
on the 1 : 18 subharmonic attractor phytoplankton and zooplankton con- 
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centrations can take at most 18 possible discrete asymptotic values (circles 
in the left panel of Fig. 4). Consequently, at sufficiently long time these in- 
teractions will lead to sharper color contrasts in the plankton distributions 
than autonomously cycling interactions, for which plankton distributions take 
continuous distributions. 

Finally, t = 60 d is far beyond the predictability horizon estimated for the 
chaotic predator-prey interactions in the right panel of Fig. 4. Correspond- 
ingly, the distributions of phytoplankton and zooplankton are nearly com- 
pletely mixed (Figs. 5-6, right panel in the bottom row), which is a result 
of the exponential divergence in time of initially nearby phytoplankton and 
zooplankton states rather than of chaotic stirring. 



4-2.3 Remarks on Other Interaction Scenarios 

We close this section briefly discussing (but not showing) the results for the two 
relevant intermediate nonautonomous biological interactions scenarios men- 
tioned earlier. These are nonautonomous interactions that take place on a 
quasiperiodic cycle and nonautonomous interactions lying on multiple periodic 
or quasiperiodic attractors. The quasiperiodically cycling interaction case is 
not different than autonomously cycling interactions as these are insensitive to 
initial conditions. The multiple attractor case, in turn, is similar to the chaotic 
interaction case at long time because the regular (periodic or quasiperiodic) 
attracting set reached by phytoplankton and zooplankton densities is sensitive 
to their initial values. 



5 Discussion 



Suppose that random velocity field fluctuations of the form u' = ^/2D(r] l (t), 

r] 2 (t)) are superimposed on u = (dyijj, — d x ip), e.g., with ip defined by (2.1). 
Here D is a constant and {t^} are Gaussian- white-noise random variables 
satisfying (%) = and (r]i('t) r 1j{ T )) = b~ijb~(t — t), where the angle brack- 
ets denote ensemble average. The expected distribution of a passive tracer, 
(9), satisfies [cf., e.g., |38| the Fokker-Plank (advection-diffusion) equation 
d t (9) + [(d) = DA (9) with Fickian diffusivity constant D, where A is 
the Laplacian operator in R 2 . Because u' is a distribution rather than a reg- 
ular function, the resulting time-dependent vector field, u + u', is no longer 
Hamiltonian — at least in a strict sense. Accordingly, KAM tori become per- 
meable to transport and the resonance islands turn leaky. In the long run, 
the mixed phase-space structures erode completely, and consequently also the 
spatial patterns in tracer distributions. 
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Standard ecological theories of pattern formation [cf., e.g.,|39| would approach 
the plankton distribution problem by replacing (4.1) with diffusion equations 
coupled by predator-prey dynamics, which is equivalent to considering (4.1) 
with a Gaussian-white-noise random velocity field. Clearly, this approach can- 
not result in a correct description of plankton advection by an oscillatory ve- 
locity field such as the one considered here because the random velocity field 
cannot produce Lagrangian mixed phase-space dynamics. Remarkably, it has 
not been until very recently £17j that standard ecological theories have been ex- 
plicitly recognized as being incapable of providing an appropriate description 
of the formation of spatial patterns in the distribution of plankton advected 
by ocean currents. In 0| advection equations like (4.1) are considered where 
the biological coupling is similar to the one considered here, except that the 
carrying capacity is set to continually relax toward a spatially nonuniform 
prescribed function (parameters are chosen so that the biological interactions 
always contain an asymptotically stable equilibrium). The velocity field in 
is a kinematic model for mesoscale turbulence, which is envisioned clS db SU. - 
perposition of randomly distributed Gaussian eddies; generation of plankton 
patchiness is mainly attributed to stirring by this velocity field. 



The traditional approach to plankton advection by turbulent ocean currents 
considers, instead of (4.1), a set of biologically-coupled advection-diffusion 
equations with eddy diffusivities chosen several orders of magnitude larger 
than molecular diffusivities. Clearly, this choice leads to a complete erosion 
of any mixed phase-space signature underlying the advection field at a much 
faster rate than molecular diffusion. More sophisticated, higher-order Marko- 
vian stochastic models have been proposed |4CJ to account for the observed 
anomalous (i.e., non-Brownian or non-random- walk) dispersion of passive trac- 
ers advected by ocean currents, which we believe is more naturally explained 
by Lagrangian mixed phase-space dynamics. Higher-order Markov models do 



not 41] give a satisfactory description of the anomalous dispersion behav- 
ior typically observed in the ocean. Other stochastic models based on non- 
Gaussian (Levy) statistics have been proposed 42||, which result in diffusion- 
type equations involving fractional partial derivatives. These so-called strange 
kinetic theories, unlike the higher-order Markovian models, have the anoma- 
lous dispersion behavior built in, but their practical utility is still not clear. 



The present work provides a basis for an alternative approach to the problem 
of pattern formation in plankton distributions in the presence of noisy ad- 
vection fields (i.e., large-scale advection fields with small-scale environmental 
perturbations superimposed) or in turbulent open ocean environments (e.g., 
mesoscale turbulence). The advection field in those environments may be de- 
composed into some mean field plus noisy perturbations (not necessarily small) 
due to environmental variability. The traditional approach would treat the 
latter as stochastic. An approach in line with the treatment given in this 
paper would consist in representing the noisy perturbation field in Fourier 
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space. Namely, through the superposition of many harmonics (e.g., Rossby- 
wave-like modes in the case of mesoscale turbulence) with random phases 
and amplitudes selected from some modeled or observed wavenumber spec- 
trum (e.g., plausibly following a power law in mesoscale turbulence). For a 
given random phase realization, the result is a velocity field with a stream- 
function that in practice is of the time-quasiperiodic form (2.3b). This allows 
for a quasi-deterministic description wherein Lagrangian mixed phase-space 
dynamics plays a central role. 



6 Summary and Conclusions 

In this paper we have investigated the effect Lagrangian dynamics on the for- 
mation of spatial patterns in the distribution of plankton near the surface 
of the ocean. Advection is assumed to be provided by a divergence-free two- 
dimensional oscillatory velocity field. Unlike the standard chaotic advection 
paradigm, in which case the oscillation of the velocity field is periodic, the 
oscillation is here assumed to be quasiper iodic, i.e., in the form of a superpo- 
sition of periodic functions with incommensurate frequencies. We analyzed in 
some detail the two-frequency case, which is suitable for describing tidal mo- 
tion in shallow seas dominated by a predominantly semidiurnal regime includ- 
ing the spring-to-neap cycle that results from superposing the semilunar and 
semisolar tidal constituents. The results, however, apply also when many more 
frequencies are considered, as would be required to model more complicated 
(turbulent) advection fields. The evolution of phytoplankton and zooplankton 
is assumed to be controlled by a set of Liouville (advection) equations coupled 
by deterministic, spatially homogeneous predator-prey dynamics of the stan- 
dard Rosenzweig-MacArthur class. In addition to autonomous interactions, 
we considered nonautonomous interactions by letting the food availability to 
prey cycle diurnally. 

The Lagrangian motion is described by a nonautonomous Hamiltonian sys- 
tem with one degree of freedom. The phase space of this motion, which is 
the physical space for advection, exhibits a mixed topological structure. This 
includes regions of chaotic motion, possibly separated by KAM tori (transport 
barriers), and populated by resonance islands (fluid particle traps). The distri- 
bution of a passive tracer is thus seen to include vigorously stirred regions and 
less stirred patches in the form of filaments and spiral-like vortices identifiable 
with the features that comprise the mixed phase space. These spatial patterns 
are subject to a quasiperiodic vibration with basic frequencies given by the 
frequencies of the time-quasiperiodic velocity field. In the particular 
lyzed here, the vibration is in synchrony with the assumed quasiperiodic tidal 
ebb and flood surpluses. This contrasts with standard chaotic advection in 
which case such a vibration is periodic. 
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The extent to which phytoplankton and zooplankton distributions are deter- 
mined by the Lagrangian mixed phase-space dynamics was shown to depend 
on the type of predator-prey interactions taking place. We considered in de- 
tail four possible interaction types: three regular (autonomous stable-spiraling, 
autonomous cycling, and nonautonomous lying on a periodic attractor) and 
one complex (nonautonomous taking place on a strange attractor). For au- 
tonomously cycling interactions and nonautonomous interactions lying on a 
periodic attractor, Lagrangian mixed phase-space dynamics leads to the for- 
mation of spatial patterns (patchiness) in the distribution of phytoplankton 
and zooplankton at all stages of the evolution. For autonomous interactions 
including an asymptotically stable equilibrium, phytoplankton and zooplank- 
ton distributions are determined by the mixed phase-space structures during 
the early stages of the evolution, before complete homogenization is reached. 
When the population dynamics is chaotic, the influence of Lagrangian mixed 
phase-space dynamics is limited to a time interval that is not too long com- 
pared to the predictability horizon. Two further relevant nonautonomous bi- 
ological interactions scenarios were also discussed in the paper. These are 
nonautonomous interactions taking place on a quasiperiodic cycle and lying 
on multiple periodic or quasiperiodic attractors. In terms of plankton distribu- 
tions, the former is not very different than autonomously cycling interactions, 
whereas the latter is similar to the long-term chaotic interaction case because 
the regular (periodic or quasiperiodic) attracting set reached by an orbit is 
sensitive to the initial conditions. 
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Fig. 5. Snapshots of the evolution of phytoplankton, which, while interacting with 
zooplankton through predator-prey dynamics governed by (4.3), is advected by the 
time-quasiperiodic tidal velocity field defined by streamfunction (2.1) with param- 
eters as in the right panel of Fig. 1. Predator-prey interaction parameters in the 
top and mid-top panels are, respectively, as in the left and right panels of Fig. 3. 
Parameters in the mid-bottom and bottom row panels are, respectively, as in left 
and right panels of Fig. 4. 
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Fig. 6. As in Fig. 5, but for zooplankton. 
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